Ultrasoft primitive model of polyionic solutions: structure, aggregation, and dynamics 
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We introduce an ultrasoft core model of interpenetrating polycations and polyanions with con- 
tinuous Gaussian charge distributions, to investigate polyelectrolyte aggregation in dilute and semi- 
dilute, salt-free solutions. The model is studied by a combination of approximate theories (random 
phase approximation and hypernetted chain theory) and numerical simulations. The calculated pair 
structure, thermodynamics, phase diagram and polyion dynamics of the symmetric version of the 
model (the "ultrasoft restricted primitive model" or UPRM) differ from the corresponding prop- 
erties of the widely studied "restricted primitive model" (RPM) where ions have hard cores. At 
sufficiently low temperatures and densities, oppositely charged polyions form weakly interacting, 
polarizable neutral pairs. The clustering probabilities, dielectric behavior and electrical conductiv- 
ity point to a line of sharp conductor-insulator transitions in the density-temperature plane. At 
very low temperatures the conductor-insulator transition line terminates near the top of a first or- 
der coexistence curve separating a high-density, liquid phase from a low-density, vapor phase. The 
simulation data hint at a tricritical behavior, reminiscent of that observed of the two-dimensional 
Coulomb Gas, which contrasts with the Ising criticality of its three-dimensional counterpart, the 
RPM. 

PACS numbers: 82.70.Dd, 64.75.Xc, 61.20.Gy, 66.10.Ed 



I. INTRODUCTION 

Ever since the pioneering work of Gouy [l|, Chapman 
[H and of Debye and Hiickel Q , electrostatic interactions 
are known to play a dominant role in determining the 
structure, dynamics and phase behavior of ionic liquids 
and solutions, as well as in governing the colloid sta- 
bility of polyelectrolyte solutions, complex biomolecular 
assemblies, and related soft matter systems. In the case 
of solutions and melts of microscopic cations and anions 
(microions), like Na + and CI - , a widely studied model 
system is the primitive model (PM) of oppositely charged 
hard spheres, which is now known to undergo a phase 
separation into a very dilute phase of mostly paired ions 
("Bjerrum pairs" Q) and a more concentrated solution 
of non-aggregated ions, for sufficiently strong Coulomb 
coupling [ll; for a review of simulation work, see [||. 

Moving to the mesoscopic scale, charged colloidal par- 
ticles (macroions) are stabilized in aqueous dispersions by 
the formation of electric double-layers of microscopic co- 
and counterions, leading to a screened Coulombic repul- 
sion between equally charged, "dressed" colloids; for an 
overview see Q • Despite this purely repulsive effective in- 
teraction between colloids, so-called "volume terms" , as- 
sociated with the self energy of individual electric double- 
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layers, induce a phase separation between dilute and con- 
centrated colloidal dispersions Q , which may be regarded 
as the highly asymmetric counterpart of the fluid-fluid 
transition of the PM [5, 6]. Recently the experimental 
and theoretical attention has shifted to the rich phase be- 
havior of such "colloidal electrolytes" , where the polyan- 
ions and polycations are highly charged, hard colloidal @ 
or nanometric particles 10 , [ll | , usually in the presence 
of added salt [13 - [l4{ . These mesoscopic electrolytes may 
thus be regarded as a generalization of the PM, where the 
pure Coulombic interactions are replaced by screened re- 
pulsive and attractive electrostatic (or Yukawa) forces. 

Another important class of complex ionic systems are 
synthetic or natural polyelectrolytes, i.e., solutions of 
charged polymer chains [l5j, where overall charge neu- 
trality is ensured by either microscopic or mesoscopic 
counterions. The latter case corresponds to a binary 
system of polymeric polyanions and polycations, which 
have been shown to aggregate into neutral or charged 
polyelectrolyte complexes (complex coacervation) , in the 
presence or absence of added salt fl6l - [l9| . Since polyions 
arc now flexible, worm-like charged objects, they can- 
not be reasonably modeled by charged hard spheres 
(like their colloidal counterparts), unless they collapse 
into quasi-spherical globules like certain folded proteins 
(e.g. lysozyme). Modern theoretical descriptions of poly- 
electrolyte complcxation are usually based on statistical 
field-theoretic formulations pol . l2lf . within a perturba- 
tive [22[ or a simulation 23 1 framework. 

In the present paper we introduce and investigate 
a simple model of polyanion/polycation aggregation in 
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polyelcctrolyte solutions without added salt. Swollen 
polymer coils in good solvent are known to interpenetrate 
easily. In fact the free energy penalty for two self-avoiding 
polymer coils to fully overlap, such that their centers-of- 
mass (CM) coincide, is of the order of twice the thermal 
energy ksT, independently of molecular weight j23 - l26| . 
A convenient coarse-grained representation of dilute and 
semi-dilute polymer coils in good solvent reduces the coils 
to ultrasoft, interpenetrating particles; the effective pair 
potentials between CM's of interpenetrating coils can be 
extracted from full monomer Monte Carlo (MC) simula- 
tions of the underlying microscopic model by a systematic 
inversion procedure p5l |26| . The resulting pair potential 
is well represented by a Gaussian of amplitude ~ 2k^T, 
and width of the order of the radius of gyration, R g , of the 
polymer coils (2|| [26| . Here we generalize the ultrasoft 
core representation to globally neutral binary solutions 
of polyanions and polycations in a dielectric continuum 
representing the solvent. The total charge of each polyion 
is assumed to be smeared over the volume of the coil ac- 
cording to a quenched Gaussian distribution of width R g , 
centered on the coil CM, ensuring that the total electro- 
static interaction between two coils of equal or opposite 
signs remains finite, even at full overlap. In particular the 
"Coulomb collapse" of oppositely charged ions, which is 
prevented by the presence of the hard core in the PM, 
is bypassed in the present model of oppositely charged 
polyelcctrolytcs by averaging over the spatial extension 
of the polyion charge distributions. 

In the following sections we investigate the pair struc- 
ture, thermodynamics and polyion dynamics of this "ul- 
trasoft primitive model" (UPM) by a combination of ap- 
proximation schemes borrowed from the theory of ionic 
liquids, and of Monte Carlo (MC) and Molecular Dy- 
namics (MD) simulations. Special emphasis is laid on the 
static and dynamic characterization of the polyion aggre- 
gation and complexation, which are expected to induce a 
conductor-insulator (CI) transition at low temperatures 
and concentrations and eventually to p hase separation 
as in the case of the PM [|, H [23, The complete 
phase diagram of the symmetric version of the UPM (the 
"restricted" UPM or URPM), which differs considerably 
from that of the RPM, will be presented in a subsequent 
paper. A preliminar y a ccount of parts of this work was 
published elsewhere j29| . 



II. THE MODEL 

We consider a binary system of N + polycations of 
charge = Z + e and iV_ polyanions of charge Q_ = 
Z_e (where e is the proton charge), moving in a dielec- 
tric continuum of dielectric permittivity e' (the "solvent" 
in its "primitive" representation) and confined to a vol- 
ume V. If n a = N a /V (a = +, — ) are the corresponding 
number densities, overall charge neutrality requires that: 



Z+n+ 



Z-n- = 



(1) 



while the total polyion number density is n = n + + 
Since the system under consideration is supposed to be 
a model for polyelectrolyte coils, the polyions are not 
point particles, but their charges are smeared over a vol- 
ume of the order of the cube of their radius of gyration 
{R g + or R g -), according to a Gaussian charge distribu- 
tion centered on the position of the CM of each polyion 
(1 < i < N = N + + N-). If r is the distance relative to 
that position, the normalized charge distribution in units 
of e of a polyion is assumed to be Z a p a (r) where 
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2-kgI 



3/2 



cxp[-r 2 /2a2] 



(2) 



with a a of the order of R ga and a = +, — ; its Fourier 
transform is 



p a (k) = J e ikl > a (r)dr = exp[-A:V/2] 



(3) 



The electrostatic potential <p a {r) generated by the 
above charge distribution obeys Poisson's equation (in 
csu) 



V 2 <p a (r) 



4nZ a e 

—Pair). 



(4) 



Taking Fourier transforms of both sides, and remember- 
ing Eq. (J3J, one finds 



ip a {k) 



4tt Z„e 



exp[-fcX/2] 



Inverse Fourier transformation leads to 

f a (r) = -7^ erf (r/V2a a 
e r \ 



(5) 



(0) 



where erf(x) is the error function. In the point particle 
limit, a a —> 0, (p a {r) reduces to Z a e/e'r, as expected. 
For finite er Q , (p a (r) remains finite as r — > 0. 

The pair potential between an a-polyion and a f3- 
polyion at a CM-CM distance r is given by 



v a p{r)= J (p a (r')Z p ep (\r-r'\)dr'. (7) 
Applying the convolution theorem, the Fourier transform 



'Jap(k) = Zpe(p a {k)p(3(k) = 



AnZaZpe 2 2 2 
— exp[-fc a af} ] 



(8) 

where o* p = (a 2 a + <rj)/2. 

Inverse Fourier transformation finally yields the set 
of three pair potentials between identical or opposite 
polyions 



v a p(r) = Q"® 13 erf (r/2cr Q(3 ) 



(9) 



This pair potential remains finite at full overlap, i.e., as 
r -> 



v afj (r) ~ u a p 

r— s-0 



C(r 6 ) 



(10) 
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where the overlap energies are 

QaQ/3 



U a /3 



(11) 



At large distances v a p(r) goes over to the Coulombic 
pair potential between point ions: 



v a p (r) 



QaQ/3 



r— ^oo CT 



1 



(12) 



In order to introduce reduced (dimensionless) physi- 
cal quantities, we define the following convenient length, 
energy and time scales 



length scale : a = 

energy scale : u = —u + 

time scale : r = 



_ Q+\Q-\ 

T 2 \ 1/2 
'+- 



(13a) 
(13b) 

(13c) 

where m is the smaller of the polyanion and polycation 
masses to_ and m+. Reduced distances, densities, times 
and temperatures are defined as 

r*=r/a < = n Q cr 3 t* = t/r T* = k B T/u 

and in the following only reduced units will be used and 
asterisks are dropped for convenience. Another, related 
dimensionless parameter is the reduced Bjerrum length 
for monovalent ions, ^b, which is identical to the Coulomb 
coupling parameter T commonly used in the characteri- 
zation of strongly coupled plasmas (see e.g. (30j ) 



T = h 



(14) 



r is directly related to the reduced inverse temperature 
/3 for a system of polyions by 



r = 



A 1 

\Z+Z\T 



(15) 



Two relevant special cases of the UPM are the follow- 
ing: 

(a) The restricted model (URPM) is the symmetric 
version where polyions arc of the same size (cr + = 
er_ = <r) and of opposite charge {Q+ = —Q- = Q)- 
This is the model introduced in our preliminary 
communication [2j|, and most of the results in the 
subsequent sections will be for the URPM. Note 
that the set of reduced units introduced above is 
different from the one we employed in Ref. (29j . 

(b) A fully asymmetric version of the model (UAPM) 
consists of polyanions (Z = \ZJ\ 1) and micro- 
scopic cations (Z + — 1), with cr_ ^> <r + . This ver- 
sion provides a simple representation of an anionic 
polyclcctrolytc in good solvent and its countcrions. 



At sufficiently low temperature and density, polyions 
and polycations of the UPM are expected to cluster into 
neutral or charged aggregates, as in the case of the PM 
[111 El SI- In the T ->■ limit, the ground state of the 
symmetric version (the URPM) is achieved by associating 
the N polyions into N/2 neutral, non-interacting pairs of 
total energy 

N 

U = --u. (16) 

This energy is finite and extensive, so that the URPM is 
thermodynamically stable according to Ruelle's stability 
criterion [33| . 

Similar considerations are expected to apply to asym- 
metric versions of the UPM, but the ground state analy- 
sis is much less straightforward, except when \Z-\ is an 
integer multiple of Z + (or conversely). To simplify nota- 
tions we consider the case where Z + = 1 and \Z- \ = Z, 
corresponding to the UAPM. One polyanion and Z coun- 
tcrions may then collapse into a neutral point cluster and 
the corresponding ground-state energy of N + /Z neutral 
clusters would be 



U a = 



N+ 
Z 



Zv + -{r = 0) + 



z(z-r 



v ++ (r = Q) 



-N + u 



1 - 



Z-l 1 



Z 2 3 / 2 



(1 + S 2 



,V2 



(17) 



where E = o"_/cr + ^> 1. 

For the ground state to be stable, Uo must be nega- 
tive. Assuming Z to scale like E 3 (volume charge), Uq is 
negative provided Z satisfies the inequality 



Z 1 > 



(z 



1 + z 2/3 



which requires Z < 22. For larger anion to cation charge 
ratios the neutral clusters are unstable, and the ground 
state can no longer be a system of neutral point clusters. 

Before presenting the results of our calculations for the 
UPM in the following sections, two remarks are in or- 
der. First, the UPM is a purely Coulombic system and 
no other force fields are involved. In particular, we do 
not include the ultrasoft entropic repulsion between inter- 
penetrating coils mentioned in the Introduction. Apart 
from obvious reasons of simplicity, this is physically jus- 
tified at low temperatures (where most of the interest- 
ing physics will be shown to occur), where Coulombic 
interactions dominate the effective entropic interactions, 
which scale like T. Secondly, a model somewhat sim- 
ilar to the UPM was used previously to investigate a 
very different Coulombic system namely a semi-classical 
Hydrogen plasma under astrophysical conditions of high 
temperatures and densities [34| . 



III. PAIR STRUCTURE AND 
THERMODYNAMICS 

The local pair structure of the UPM is characterized 
by the three partial pair distribution functions g++{r), 
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(?H (r), and g (r); the corresponding pair correlation 

functions h a p (r) = g a p (r) — 1 go to zero at long distances 
in the disordered fluid phases. Their Fourier transforms 
h a p(k) are related to those of the direct correlation func- 
tions c a p(k) by the set of three coupled Ornstein-Zernike 
(OZ) relations [IH; solving the latter for the h a p(k), and 
introducing the concentrations x a = n a /n, we find 



i++(k) 



1 



D{k) 

h + -(k) = h-+(h) = 
1 



c ++ (fc) - x-A(k)) 
c+-(k) 



■(*) 



D(k) 



D(k) 
c—(k)-x + A(k)} 



(18a) 
(18b) 
(18c) 



where 



A(fc) = c ++ (fc)c__(fc) 



-(*) 



D(k) = [1 - x+c ++ {k)] [1 - x_c__(fc)] 

— X+X-C^ (k) 



(19) 



and all Fourier transforms are dimcnsionless, i.e., f(k) = 
n J exp(ikr)/(r)dr. 

The partial structure factors S a p(k), which are mea- 
surable by X-ray or neutron diffraction experiments, are 
directly related to the h a p(k) through (3f| 



S a p{k) = — (/3kP_k) = x a$a/3 + X a Xph a p{k) (20) 

where the Fourier components of the local density oper- 
ators are 



k = ^ expjik ■ v ia ). 



(21) 



It is convenient to introduce the linear combinations 
corresponding to total number (N) and charge (C) den- 
sities 



N + i - 
Pk =Pk+Pk 

Pk = z +p£ + z -pk 

and the related structure factors 

SNN(k) = ^(p£p N k ) = £E« fc ) 

a p 

a p 

scc(k) = ^{ P <ip%) =Y,Y. z « z M k ) 



(22a) 
(22b) 



(23a) 
(23b) 
(23c) 



a p 



The charge-charge structure factor obeys the Stillingcr- 
Lovett limit, valid for a conducting medium 31 1 



k^o k 2 Z 2 



(24) 



where Z 2 — x+Z 2 ; + x_Z 2 _ and kd is the inverse Debyc 
screening length 



_2 



V D- 



Aim+Zlln + Aim-Z 2 l 



-LB- 



(25) 



Knowledge of the pair structure gives access to a num- 
ber of thermodynamic properties via standard relations 
[H| ■ Using the dimensionless variables defined in Section 
HTl the reduced excess internal energy per polyion is given 
by 

..ex W«* 



N 



= 2-nnT / {x\Z%h ++ {r)exi{r/2a + ) 
Jo 

+2x + x_Z + Z_h + _{r)ei{{r/2<T + _) 
+x 2 _Zlh—{r)cri(r/2<7-)}rdr (26a) 



I" Jo 1 



x%Z 2 + h ++ (k)e 



+2x+x-Z+Z-h^ (fc)e 



-K (J i 



+x z _Z±h__(k)e 



dk 



(26b) 



where the transition from the first to the second of the 
above equations has been achieved by using Parseval's 
theorem. The wavenumbers k in Eq. (|26b[) arc dimen- 
sionless, i.e., expressed in units of cr _1 = cr^i. 

Similarly, the dimensionless equation of state (with P 
the pressure) is given by: 



BP u cx 

= 1 + ^r 

n 6 

2r0F 



• T + z + 
0+ 



h ++ (r)e 



-r 2 /<lo 



h+_(r)e- r2/4 <- 



2x+X-Z+Z- 
+ — ± — 



(27a) 



-f £ { X \Z\o\h ++ (k) C - k2 < 

+2x + X-Z + Z^al_h + -(k)e- k2,7 +- 

-(fc)e- fc2<T -}fc 2 dfc (27b) 



-x 2 Z 2 o 2 _h-( 



The dimensionless isothermal compressibility xt is deter- 
mined by the small fc-limit of the number-number struc- 
ture factor: 

nk B T XT = (^) = lim S NN (k). (28) 
\ on J T fc-H) 

The entropy S, Helmholtz free energy F or chemical 
potentials p + and /x_ cannot be expressed in terms of the 
pair distribution functions alone. The dimensionless ex- 
cess Helmholtz free energy per polyion, / cx = /3F cx /iV, 
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along an isochore can be obtained by standard thermo- 
dynamic integration, starting from the high temperature 
(/3 = 0) limit, where / cx = 0, according to 



f cx W,n) 



dff 
«° x (/?»f- 



(29) 



Equivalently, the chemical potentials of the polyions may 
be calculated by the Kirkwood charging process [36[ 
whereby a test polycation or polyanion is gradually cou- 
pled to the N polyions of the system by varying a cou- 
pling parameter A from (non interacting test particle) 
to 1 (fully interacting test particle): 



% j dX 
Jo 



^ xpv a p (r)h af} (r; X)dv 





(30) 



where h a p(r;X) is the pair correlation function between 
the test particle of species a and the bath particles of 
species /3 corresponding to a pair potential Xv a p(r). 

Because of the absence of strong short-range interac- 
tions and the long-range nature of the Coulombic interac- 
tions, a "natural" approximation for the pair correlation 
functions is provided by the random phase approxima- 
tion (RPA), which is expected to be very accurate at 
high temperatures (T > 1) and high densities. The RPA 
amounts to setting the direct correlation functions equal 
to their asymptotic limit, for all inter-particle distances 
r, namely 



r 

c a p(r) = ~/3v al3 (r) = -Z a Z l3 -cd(r/2a al3 ) 

r 

~ />\ n~ r,\ ATrZ a ZRTn _i.3 a 
c a p{k) = -pv a0 {k) = ^ — e k a ^ 



(31a) 



(31b) 



Substitution of Eq. ()31b|) into Eqs. (g5a| to (JT5J) leads 
to simple analytic expressions for the h a p{k), while 
Eqs. (|23ap to (|23c[) yield the following explicit expres- 
sions for the structure factors 



It is easily verified that lim/^o SMN{k) — 1, so that, 
according to Eq. (|2"8)l the isothermal compressibility is 
that of an ideal gas, a well-known deficiency of the RPA. 
Similarly, lim^o SNc{k) = 0, so that the RPA predicts 
that charge and number density fluctuations are decou- 
pled in the long wavelength limit. After some algebra the 
charge-charge structure factor (|32c[) may be cast in the 
convenient form 



Scc{k) 



Z 2 k 2 + x+x- 


2 i K £j 


' e -k 2 a\/2 _ e -k 2 a 2 _/2 


Z 2 k 2 + k 2 d 


x+Z\c- k 


2 <+ X -Z 2 c- k2 °-~ 





Z 2 



(34) 

which clearly satisfies the Stillinger-Lovett condition 

In the remainder of this paper we focus on the symmet- 
ric (restricted) version of the UPM, namely the URPM 



Z: n A 



7i- = n/2; (j-i 



er). In this case, v++(r) = u__(r) = —V-\ — (r), and hence 

<7++(r) = g (r), so that there are only two independent 

pair distribution functions and structure factors. Fur- 
thermore, all the relations above simplify considerably 
for the URPM. For instance, S NN {k) = 1, S NC {k) = 
and Scc{k) reduces to 



Scc(k) 
Z 2 



k 2 



K D e_fe2c 



(35) 



For point ions, Scc(k) goes over to the standard Debye- 
Hiickcl form. Note that, within RPA, h ++ (r) = 
— h + _(r) = h cc {r)j2. RPA estimates of the ther- 
modynamic properties are obtained by substituting the 
above expressions for the pair correlation functions in 
Eqs. ([26b| and (|27b| . The exact Kirkwood expression 
pop for the excess chemical potentials /i+ = /i_ = fi 
reduces to 



P» e *=2j dX I Mr)hcc(r;\)dr 



1 1 



2 (2tt) 



d\ / /3v(k)h C c{k;X)dk (36) 



'iVJV 



(k) 



1 



SNc(k) 



k 2 D{k) 
-\-2x+X-Z+Z-C 



x\Z\&~ k a + +x 2 _Z 2 _e- k ~ 

-k 2 o\_ 



x 2 + Z%e~ 



Scc{k) = Z 



k 2 D(k) 
+x + X-Z+Z-(Z+ + ZJ)<T k 



+ x 2 _Z 3 _e~ 



2 2 
+ - 



k 2 D{k) 



x 2 + Z 4 + c- k 



+2x + x-ZiZ 2 _e- k 



where according to (|19p and (|31b|) 



D(k) = 1 



Z 2 k 2 



x+Z 2 + c- k ' (T + + X -Z'ie 



2 n -k* 



(32a) 



(32b) 



+ +xlZ A _e- k *°- 



(32c) 



(33) 



and 



where v(r) = v ++ (r), /3v(k) = ""^ - e 
hcc{r) = h ++ (r) - h+-(r). 

Using the coupled OZ relations for the correlation func- 
tions of a test particle partially coupled to the bath par- 
ticles one easily shows that 



(37) 

-2nf3v(k) and 

c cc {k]X) = -2Xpv(k) = Xccc(k) within RPA. Hence 
hcc(k; X) = Xhcc(k) so that the integration over A in 
Eq. (|3l)l) is trivial leading to 



h C c(k; A) = 
where c C c{k) = c ++ (fc) 



ccc(k; A) 
1 - \c CC {k) 

c+-(k) = 



W = 7 



4 (2tt) 3 
TZ 2 p °° 



2n 



h C c{k)pv{k)dk 
hcc{k)e~ k2,j2 dk = u^. (38) 
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The last equality follows from Eq. (|26b[) adapted to the 
symmetric case. The result that the excess chemical po- 
tential is equal to the excess internal energy per ion is of 
course only true within RPA. 

It is easily verified from (|38|) and (|35|) that the low 
temperature limit of the RPA internal energy per ion 
coincides with the ground state energy (|16|) 



lim 

T-s-0 



u 



RPA 



(T) = _N = Uo 
N 2 U N ' 



(39) 



In the high temperature limit the leading contribution to 
the reduced excess energy per ion reduces to the Debye- 
Hiickcl limiting law for point ions 



hm — 

T->oo N 



inn 



(40) 



In the zero temperature limit the pressure is negative 

pRPA 



lim 

T-yO 



u 
12' 



(41) 



Among the standard integral equations for the pair 
structure the hypernetted chain (HNC) equation is 
known to be well adapted to Coulombic fluids [35[ . HNC 
theory supplements the coupled OZ relations linking the 
total and direct correlation functions h a p{r) and c a p{r), 
by the closure relations 

h a p{r) + 1 = g a p(r) = cxp{-/3u Q/3 (r) + h a/3 (r) - c Q/3 (r)} 
= cxp{Va(r) - Ac Q/3 (r)} (42) 

where Hs.c a p{r) = c a p(r) + fiv a p{r) are the short-range 
parts of the direct correlation functions, which are ex- 
pected to vanish rapidly for large r (within the RPA, 
Ac a p(r) = 0). In the symmetric case under consider- 
ation here, /i ++ (r) = h (r), c ++ (r) = c (r), and 

the three coupled OZ relations reduce to two decoupled 
equations for the number-number and charge-charge cor- 
relation functions hjyN{r) = h ++ (r) + /i + _(r), hcc( r ) = 
/i ++ (r) — h-\ (r); in fc-space the two relations read 



h aa (k) — 



c aa (/c) 



1 



«(*)' 



N, C. 



(43) 



HNC theory has the advantage that the excess chemi- 
cal potential can be calculated from a knowledge of the 
pair correlation functions alone [35[ . The set of Eq. (|4"2")l 
and (|43[) are solved numerically by an iterative Picard 
method, taking particular care of the long-range tails of 
the direct correlation functions in r- and fc-space. Un- 
fortunately, the range of thermodynamic conditions for 
which the numerical procedure converges is limited to in- 
creasingly higher temperatures as the density is reduced. 
The applicability of the HNC closure is therefore limited 
to relatively high density and temperatures (see below). 

We have compared the RPA and HNC predictions for 
the pair structure of the URPM to MC results along sev- 
eral isochores. The details of our simulations are de- 
scribed in Appendix |XJ In the following we focus on 
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FIG. 1. Total correlation functions h++(r), h+-(r) and 
hcc{r) from MC simulations (symbols) and RPA (full lines) 
for two different state points: (a) n = 0.0035, T — 0.63 and 
(b) n = 0.35, T = 0.25. 
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FIG. 2. Pair distribution functions along the isochore n = 
0.0035 for three different temperatures (as specified) from MC 
simulations (symbols) and RPA (full lines): (a) gj, — (r) and 
(b) s++(r). 



the two isochores n = 0.0035 and n = 0.35, which are 
representative of the system's behavior at low and high 
density, respectively. Note that we do not investigate 
here the regime of very low densities, in which special 
simulation techniques must be employed to ensure er- 
godicity (37} • Data along the isochores n — 0.0035 and 
n = 0.35 are shown in Fig. Q] to 21 respectively. Fig- 

ure[T] shows h ++ (r) 1 h-\ (r) and hcc( r ) for n = 0.0035, 

T = 0.63 and n = 0.35, T = 0.25. Since these tempera- 
tures are below the no-solution line of HNC, we only re- 
port the predictions of RPA. At the higher density, RPA 
is seen to be very accurate for hec ( r ) ■ although the RPA 
symmetry (r) = — /i55 (r) is broken in the Simula- 
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FIG. 3. Pair distribution functions along the isochore n = 
0.35 for three different temperatures (as specified) from MC 
simulations (symbols), RPA (full lines), and HNC (dashed 
lines): (a) g+-(r) and (b) g++(r). 



FIG. 4. Structure factors along the isochore n = 0.0035 for 
three different temperatures (as specified) from MC simula- 
tions (symbols) and RPA (full lines): (a) Scc(k) and (b) 
5jvjv(fe). 



tion data. The discrepancies between RPA and simula- 
tion data are only slightly more pronounced at the lower 
density (n — 0.0035, T = 0.63), but they increase rapidly 
as T is lowered due to strong Coulomb correlations (see 
below). 

A comparison between theoretical predictions and MC 
data at even lower temperatures is made in Fig. [5] along 
the isochore n = 0.0035, and along the isochore n = 0.35 
in Fig. [3] Strong pairing is evident from inspection of 
the data for (r) at the lower densities and temper- 
atures as expected. Indirect evidence of pairing is also 
provided by the observation that at the lowest temper- 
atures, along the isochore n = 0.0035, <?++(r) > 1 as 
r — > despite the Coulomb repulsion between equal sign 
polyions; this apparent attraction can be rationalized by 
the formation of tight anion/cation pairs: the anion of 
a given pair will be preferentially situated between its 
cationic partner and the cation of another pair, thus fa- 
voring the cation-cation approach and the possible for- 
mation of trimers and higher order clusters. RPA is un- 
able to properly account for strong ion pairing, both at 
high and low density. The agreement of HNC predic- 
tions with the MC data at high density (see Fig. [3J is 
good at high temperature but is seen to deteriorate as 
the temperature is lowered towards the no-solution point 
(corresponding to T = 0.044 at n = 0.35). Also note 
the unphysical predictions of RPA (g++(r) < 0) at low 
density and temperature. 

Figure HJ-a shows a comparison between the RPA pre- 
dictions and MC data for the charge-charge structure fac- 
tor Scc{k) along the low density isochore n = 0.0035. 
While the agreement is very good at the higher temper- 
atures (T > 0.2), it deteriorates rapidly as T is lowered. 
At the lowest temperatures (T < 0.05) the MC data do 
not satisfy the Stillinger-Lovett condition (|2~4"|) . while the 



RPA result j34| obeys the condition by construction. As 
shown in |29j . a least squares parabolic fit to the low-fc 
MC data [S C c{k)/Z 2 ~ k 2 /n 2 } yields an effective in- 
verse Debye screening length k which is systematically 
larger than ftrj [defined in Eq. (|2~5j) ]. The violation of the 
Stillinger-Lovett limit suggests that the URPM no longer 
behaves as a purely ionic fluid at low temperatures, due 
to ion pairing; we will return to this important issue in 
Sec. IIV| where clustering and dielectric response will be 
analyzed. A similar analysis of Scc(k) along the high 
density isochore (see Fig. [5]) shows excellent agreement 
between RPA and MC data down to T = 0.003, with the 
effective K consistently close to k^, confirming that the 
URPM remains ionic throughout: at high densities an- 
ion/cation pairs are short-lived due to the frequent over- 
lap of neighboring pairs. 

The results for SNN{k) are shown in Fig. HJ-b and 
Fig. [5]4d. The density-density structure factor S^Nik) 
is seen to increase rapidly as k — > along the isochore 
n = 0.0035, signaling the proximity of a spinodal line as 
T is lowered. At the spinodal line associated with a phase 
separation, SNN(k) diverges in the long wavelength limit. 
Note that, within RPA, SataK^) = 1. The predictions of 
HNC along the high density isochore indicate a mild in- 
crease at small k, consistent with the simulation data. 
However, the integral equation approach admits no solu- 
tion in the vicinity of phase coexistence — a well-known 
defect of the theory [38l . l39j . All these findings clearly 
indicate that more sophisticated theoretical tools should 
be employed to study the low temperature, low density 
regime of the model. 

The RPA, HNC and MC data for the pair struc- 
ture may be used to compute the excess internal energy 
and the equation of state of the URPM via Eqs. (|26al) . 
()26b|) and (|27a[) . (|27b|) . and the isothermal compress- 
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FIG. 5. Structure factors along the isochore n = 0.35 for 
three different temperatures (as specified) from MC simula- 
tions (symbols), RPA (full lines) and HNC (dashed lines): (a) 
Scc(k) and (b) S NN {k). 




FIG. 6. Excess internal energy per particle U ex /N as a func- 
tion of temperature from MC simulations (symbols) and RPA 
(full lines) along the isochores n = 0.0035 (squares) and 
n = 0.35 (circles). 



ibility via Eq. (|28l) . The excess chemical potential fol- 
lows from Eq. ([55]) . leading to the explicit result 
within the RPA. Representative comparisons are made 
for u cx (n, T) = U cx /N in Fig. [5] along the isochores 
n = 0.0035 and n = 0.35, and for the total, dimensionless 
chemical potential /3/j, in Fig. [7J Simulation results for 
the density dependence of (3fi have been obtained from 
grand-canonical MC simulations using biased pair inser- 
tions [4p| . Also included are results for /3/i obtained from 
standard NVT MC simulations using the Widom inser- 
tion method [4l[ . The agreement between RPA and sim- 
ulation results is good only at sufficiently high temper- 
ature and density. We note that for these temperatures 
and densities, RPA- and HNC-data agree within symbol 
size. For completeness, we note that the Hclmholtz ex- 
cess free energy per ion can finally be obtained from the 



FIG. 7. Dimensionless chemical potential /3fi as a function of 
the density n along several isotherms (as labeled) obtained 
from the Widom insertion method in NVT-MC simulations 
(open circles), grand canonical MC simulations (full circles) 
and RPA (full lines). 



thermodynamic relation 

f3F™ 



r 



N 



(44) 



IV. 



CLUSTERING AND DIELECTRIC 
RESPONSE 



The structural data and thermodynamic properties re- 
ported in section HTT1 point to strong clustering of anions 
and cations at low temperatures and densities. Pair- 
ing of oppositely charged ions, as well as the formation 
of larg er aggregates are well documented for the RPM 
[3, |27| . [HI, |33, , where clusters carry large multipole 
moments (dipoles in the case of pairs) due to the ionic 
cores which induce charge separation. In the case of the 
URPM opposite charges can overlap, so that clusters do 
not carry permanent multipoles, but are polarizable en- 
tities. Their mutual interactions are thus expected to be 
much weaker than in the case of the RPM. 

Clusters can be properly defined only at sufficiently 
low densities, and even then there is always some arbi- 
trariness in their definition. We have used a standard 
geometric definition of m-clusters, namely that m ions 
form an m-mer if each ion lies within a given distance 
r c of at least one other ion in the cluster. The cut-off 
r c is taken to be typically close to 1.0 (in units of a), 
corresponding to a situation where the charge distribu- 
tions of polyions touch (cf. Eq. ([5])), and we have exam- 
ined the sensitivity of the cluster analysis to variations 
of r c . We have thus determined the cluster distribution 
functions, i.e., the fractions P(m) of m-mcrs, averaged 
over all configurations generated in MC and MD simula- 
tions. Typical examples of the resulting histograms are 
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FIG. 8. Cluster distribution function P(m) for (a) T = 0.125 
and (b) T = 0.0125 along the isochore n = 0.0035. 



shown in Fig. |5] for n = 0.0035 and two temperatures. At 
the higher temperature (T = 0.125) monomers (i.e., free 
ions) are by far the majority species, and P(rn) is seen to 
decrease monotonically and rapidly as m increases. At 
the lower temperature (T = 0.0035), pairs constitute the 
most common m-mers (with a fraction P(2) very close to 
1), and the variation with m is alternating, with neutral 
m-mers (pairs, tetramers, etc.) being much more prob- 
able than charged m-mers (corresponding to odd values 
of m). 

The variation of P(m) as a function of temperature, 
for 1 < m < 4, is illustrated in Fig. |9] for three different 
choices of the cut-off r c . Although there are significant 
quantitative differences between the three sets of data, 
the trends are the same. As expected the fraction P(l) 
of monomers is vanishingly small at the lowest tempera- 
tures; it increases rapidly towards 1 for T > 0.05. Simul- 
taneously the fraction P(2) of dimers drops dramatically 
as T increases; the fraction P(3) of trimers is negligi- 
bly small throughout, while the fraction of tetramers is 
non negligible at the lowest temperatures. These trends 
agree with intuitive expectation, but our analysis pro- 
vides a quantitative estimate of the temperature range 
associated with the onset of pairing. The onset of pair- 
ing along the isochore n = 0.0035 is seen to be rapid, but 
continuous. A similar cluster analysis carried out along 
several low density isochorcs shows that the temperature 
at which P(m — 2) ~ P(m = 1) ~ 0.5 drops as the 
density n increases. 

Using MD simulations, we have estimated the mean 
cluster lifetimes r m for m-mers with 1 < m < 4 as 
functions of temperature, along the isochore n = 0.0035. 
The lifetime is defined as the minimum time span dur- 
ing which a cluster is formed by exactly the same set of 
particles, and is thus bounded from below by the typi- 
cal collision time between clusters. Results are shown in 
Fig. 1101 The monomer lifetime is seen to increase from 




FIG. 9. Fractions of selected cluster types P(m) as functions 
of temperature along the isochore n = 0.0035: P(m = 1) 
(squares), P(m — 2) (circles), P(m = 3) (triangles), and 
P(m = 4) (reversed triangles). The three panels show results 
obtained using the following values of the cut-off r c : (a) r c = 
1.4, (b) r c = 1.0, and (c) r c = 0.7. 
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FIG. 10. Cluster lifetime r m at n = 0.0035 as a function 
of temperature for selected cluster types: m = 1 (squares), 
7Ji = 2 (circles), m = 3 (triangles), and m — 4 (reversed 
triangles). 

Ti ~ 2 up to Ti ~ 15 as the temperature rises, while the 
opposite trend of the dimer lifetime is seen to be much 
faster, with T2 = 10 3 at the lowest T, dropping rapidly to 
Ti ~ 5 at T ~ 0.05, the temperature at which t\ and T2 
cross. The lifetimes T3 and T4 of trimers and tetramers 
are even smaller at high T, although they tend to increase 
in the low-temperature, paired regime. 

At the lower temperatures, where the lifetimes ti of 
pairs are long, it makes sense to consider the systems as 
made up of three "species" , namely free polyanions, poly- 
cations, and neutral polyanion-polycation pairs ("chem- 
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FIG. 11. Distribution functions g m ,n{r) between the centers 
of mass of selected cluster types for n = 0.0035, T = 0.025: 
(a) monomer-monomer <?i,i(r), (b) monomer-dimer gi,2(f), 
and (c) dimer-dimer 32,2 (r). 



ical picture"). Under those conditions, one may deter- 
mine monomer-dimer and dimer-dimer distribution func- 
tions (71,2 (?") and 32,2 ( r ) between isolated monomers and 
dimers that are not part of larger clusters, in MC or MD 
simulations. Examples arc compared to the monomer- 
monomer pair distribution function gi t i(r) in Fig. II H 
for n — 0.0035 and T — 0.025, i.e., close to the pair- 
ing transition. gi,i(r) is dominated by the strong anion- 
cation attraction, leading to a pronounced peak at r = 1. 
51,2 {f) and <72,2 ( r ) are seen to exhibit modest peaks at 
slightly larger distances pointing to rather weak correla- 
tions. Note that all three pair distribution functions van- 
ish for r < r c = 1.0 by construction, because for shorter 
CM-CM distances the two monomers would be identi- 
fied as a single dimer, a monomer and a dimers would 
be identified as a single trimer, and two dimers would 
be identified as a single tetramer, implying an apparent 
infinite barrier [3l| . 

In view of the low density, the pair distribution 
function data may be inverted to determine effective 
monomer-monomer and monomer-dimer pair potentials 
^1,2(0 and V2,2{f) by a simple Boltzmann inversion 



( (r) = -k B T ha. g n , m (r) 



(45) 



The resulting effective potentials are pictured in Figs.[T2l 
and 1131 As explained above, the short-range repul- 
sion is an artifact linked to cluster identity. Beyond 
r = r c the monomer-dimer and dimer-dimer potentials 
are seen to be weak and attractive. In view of the sta- 
tistical uncertainties, the asymptotic behaviors at large 
distances, illustrated by the insets in Figs.fT2landfT3l are 
in reasonable agreement with the theoretical predictions 
Vi2(f) ~ l/?' 4 and 1)2,2(1") ~ 1/r 6 for isolated monomer- 
dimer and dimer-dimer pairs (see Appendix B). These 



FIG. 12. Dimensionless, effective dimer-dimer potential 
Pv2,2(r) for n = 0.0035, T = 0.025. Inset: effective potentials 
on a double-logarithmic scale. The dashed line indicates the 
asymptotic behavior expected at large r (see Appendix [B}. 
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FIG. 13. Same as Fig. ll2l but for the monomer-dimer effective 
potential f3vi,2(r). 



are of course the same asymptotic behaviors as for van 
der Waals-London dispersion interactions between polar- 
izable atoms or molecules, but the prcfactors are propor- 
tional to k^T in the present, purely classical case, while 
they are independent of temperature in the atomic case, 
since they result from quantum averages over the elec- 
tronic ground state. 

Another diagnostic for clustering is provided by the di- 
electric response of the URPM. We have shown in Section 
lllll that at high temperatures and densities the system be- 
haves as a conductor, obeying the Stillinger-Lovett per- 
fect screening condition Eq. (f24)). However, as illustrated 
in Fig. SI the Stillinger-Lovett condition is violated at low 
densities and temperatures pointing to a dielectric, rather 
than conducting behavior. According to static linear re- 
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FIG. 14. Dielectric order parameter (e — l)/e (filled circles) 
obtained from Eq. (|49[) and fraction of free ions P(m = 1) 
(crosses) at n = 0.0035 as functions of temperature. Inset: 
enlarged view of the low-temperature behavior of e. 



sponsc theory 35] , the charge-charge structure factor is 
related to the /c-dependent dielectric response function 
by 



1 



1 



4 S cc (k) 
k 2 Z 2 ' 



(46) 



Taking the k — > limit, where Scc(k) ^ Z 2 k /k 2 we 
arrive at the following expression for the macroscopic di- 
electric permittivity e 



- = lim — 

e fc->o e{k) 



_ i K D 
K 2 



(47) 



In a conducting medium k = kd, and e — > oo as expected. 
At sufficiently low T and n, all polyions are paired and 
the MC data show that k > kd , so that e takes on a finite 
value characteristic of an insulating, dielectric medium 



oo > e 



1 



1 - 4/k ; 



> 1. 



(48) 



Equivalently, e may be estimated form Kirkwood's fluc- 
tuation formula [42| adapted to simulations carried out 
under metallic boundary conditions [43[ 



An (|M| 2 )-|(M)P 
3k B T V 



(49) 



where M = QiYi is the total dipole moment of the 
system. 

From the permittivities estimated using Eq. (|49[) it is 
possible to calculate a "dielectric order parameter", de- 
fined as (e— l)/e, which equals 1 in the conducting phase 
and is close to in the insulating phase. The values of 
(e — l)/e obtained from MC simulations along the iso- 
chore n = 0.0035 are shown in Fig. [TJ] together with 



the corresponding percentage of monomers P{m = 1). 
The dielectric order parameter signals a sharp transition 
around T ss 0.04, which correlates with the rapid increase 
of the fraction of free ions. In the inset of Fig. [THwe see 
that e remains nearly constant up to T ~ 0.025, beyond 
which e increases sharply towards a high temperature 
limit close to e ~ 10 3 (cf. Ref. [29]), typical of a finite 
conducting medium (for an infinite conductor, e — > oo). 
The correlation between (e— l)/e and P{m = 1) is strik- 
ing and validates the picture of a transition from a low 
temperature dielectric (insulator) state to a high tem- 
perature ionic (conductor) state driven by the break-up 
of ion pairs. The "conductor-insulator" (CI) transition 
is seen to be fast, but continuous, which may well be a 
finite size effect. Simulations on larger systems and a fi- 
nite size scaling analysis will be required to find out if the 
transition becomes discontinuous in the thermodynamic 
limit. 

The above analysis has been repeated along several low 
density isochores, and the locus of points in the (n,T) 
plane where P(m = 1) = 0.5 yields an estimate of pair- 
ing (or CI) transition line which is shown in Fig. [TSJ In 
the portion of the phase diagram studied herein, the tran- 
sition is seen to take place at lower temperatures as the 
density increases, contrary to what could be expected 
from a simple "chemical equilibrium" picture [27] and to 
what is found in the RPM (H. 



V. POLYION DYNAMICS 

Time-dependent correlation functions are readily cal- 
culated by MD simulations, and provide a quantitative 
characterization of single-particle and collective ion dy- 
namics as well as giving access to linear transport coeffi- 
cients, like ion mobility and electrical conductivity. The 
latter provides an unambiguous diagnostic of a CI tran- 
sition. We note that Brownian Dynamics simulations 
would be more appropriate for the system at hand, since 
they account for solvent effects. However, since we are 
interested in qualitative, rather than quantitative aspects 
of the dynamics, we prefer to employ here the more ef- 
ficient MD simulations, which are capable of exploring 
longer time windows. 

Let r(t) and v(t) denote the position and the velocity 
of any given polyion at time t; the normalized velocity 
auto-correlation functions of the anions and cations are 
identical for the symmetric URPM and defined by 



Z+(t) = Z_{t) = Z(t) 



<v(t) • v(0)> 



(50) 



where v(0) denotes the initial velocity, (v 2 ) = "Sk-sT/m, 
and statistical averages arc taken along the trajectories 
of individual ions, and over all N ions. Examples of MD- 
generated correlation functions Z(t), for n = 0.0035 and 
four temperatures, are shown in Fig. [151 At the higher 
temperatures, Z(t) is seen to decay essentially exponen- 
tially. However, the relaxation time of Z(t) increases 




FIG. 15. Velocity auto-correlation function Z(t), as denned in 
Eq. (150[) , at n = 0.0035 for different temperatures (as labeled). 
Inset: self diffusion constant D evaluated via Eq. (I52|l as a 
function of temperature. 



FIG. 16. Mean square displacement (|M(t) - M(0)| 2 > of the 
total electric dipole M(£) as a function of time at n = 0.0035 
for several temperatures (as labeled). The dashed line in- 
dicates the asymptotic behavior expected for a conducting 
system. 



with increasing T . This counter-intuitive result can be 
explained in terms of the ultrasoft nature of the poten- 
tial: at higher T collisions become less and less effective 
in decorrclating the ion velocities, since particles barely 
feel each other's influence. Thus, the "effective" collision 
time increases with increasing T. The relaxation changes 
dramatically at the lower temperatures, where Z(t) de- 
cays much more slowly after an initially strong oscilla- 
tory regime. This striking behavior may be rationalized 
in terms of anion-cation pair formation. The oscillations 
correspond to the vibrations of the ions relative to the 
CM of a "bound" pair. The measured reduced angular 
frequency, u> ~ 0.6 is close to the frequency of the rel- 
ative motion of two oppositely charged ions within the 
harmonic potential (see Eq. (fTD))) 



u+_(r) 



12a 2 



(51) 



i.e., w_| = yi/3 — 0.58, independent of temperature. 

The slow decay of Z(t) at long times may be associated 
with the motion of the long-lived pair within which the 
ion is bound. Due to the weakness of the interaction 
between pairs (cf. Fig. [12]) , the motion of neutral pairs 
is nearly ballistic, resulting in the very slow relaxation 
shown in Fig. 1151 The self-diffusion constant of the ions, 
D + = D- = D may be calculated b y in tegrating the 



velocity auto-correlation function (|50|) [35|J, or, more ac- 
curately, from the asymptotic slope of the mean square 
displacement of an ion from its initial position, according 
to Einstein's relation wfj 



D = lim 

t— >OC 



(|r(f)-r(0)p 
6t 



(52) 



The variation of D with temperature is shown in the in- 
set of Fig. [T21 D is seen to first drop as T decreases, which 
is the usual behavior observed in "normal" liquids before 
increasing sharply for T < 0.08 and go through a pro- 
nounced maximum around T ~ 0.03 below which D de- 
creases again towards zero. This unusual non-monotonic 
behavior is obviously a direct consequence of pairing: as 
argued earlier pairing leads to a system of nearly non- 
interacting entities which move almost freely, thus ex- 
plaining the sharp rise in mobility. Once pairing is com- 
plete and no free ions are left, the diffusion constant will 
drop again with temperature. 

We now turn to the electrical conductivity a e per unit 
volume of the URPM. a e is determined by the asymptotic 
slope of the mean square displacement of the total electric 
dipole M, according to the generalized Einstein relation 



1 



Vk B T t- 



lim 



\M(t) - M( 
6i 



(53) 



The corresponding ion mobility is [i = D/k^T. 



which is equivalent to the familiar Green-Kubo relation 
linking a e to the electric current auto-correlation function 
[35l ]; in view of the large statistical uncertainties of the 
latter, as generated in MD simulations, the Einstein rela- 
tion (|5"3"|) is preferable to estimate <r e . Data for the diffu- 
sion c(t) = (|M(t) - M(0)| 2 ) of the total dipole along the 
isochorc n = 0.0035 arc shown in Fig. [TCJ Note that ions 
must be allowed to leak out of the periodically repeated 
simulations cell when calculating M(i), to allow for dif- 
fusion in an unbounded volume. The long-time slopes of 
the c(i)-curves yield <r e . Figure (THl shows that the asymp- 
totic linear regime is rapidly reached at the higher tem- 
peratures, while at the lower temperatures the regime is 
only reached after reduced times t > 10 3 . At the lowest 
temperature investigated, c(f ) appears to tend to a con- 
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FIG. 17. Electrical conductivity a e , evaluated via Eq. (|53[) . 
as a function of temperature along different isochores (as la- 
beled). T c indicates the estimated temperature where a e van- 
ishes, according to power-law fits a e ~ (T — T c ) v with u — 1.2 
(full lines). The estimated uncertainty on v is ±0.02. 



stant, i.e., its slope and hence <r e , are zero, correspond- 
ing to an insulating state. Note that c(t) exhibits os- 
cillations at the lower temperatures reminiscent of those 
observed in Z(t) (cf. Fig. I15| . of comparable frequency, 
which are once again ascribable to pairing. We have re- 
peated the electrical conductivity analysis along several 
low density isochores. The data for the temperature de- 
pendence of <7 e are shown in Fig. |T7l The conductivity 
is seen to drop to zero at increasingly low temperatures 
as the density increases. Assuming a power-law depen- 
dence a e ~ (T — T c y with a fixed value of v ~ 1.2, least 
squares fits provide a set of "critical" temperatures T c at 
which <r e (T) appear to vanish (indicated by the labels in 
Fig. [T7|) . These results provide an estimate of a CI tran- 
sition line in the (n, T) plane. This line is compared in 
Fig. [18] to the pairing transition line which provides an- 
other estimate of the location of the CI transition. The 
two lines are seen to be roughly parallel and reasonably 
close. We also include in Fig. [18] the pairing transition 
line reported in Ref. [2!|, in which different values of 
r c and of the "critical" fraction of free ions P(m = 1) 
were used. As suggested by Fig. [21 different choices of 
P(m = 1) would correspond to different values of a "crit- 
ical" dielectric permittivity e. 



VI. DISCUSSION AND CONCLUSIONS 

We have introduced a model of interpenetrating 
polyanions and polycations, carrying extended, continu- 
ous charge distributions. We expect this ultrasoft primi- 
tive model (UPM) to be relevant for the study of the ag- 
gregation of oppositely charged polyelectrolytes in good 
solvent. Most of the results presented in this paper 
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FIG. 18. Current estimate of the phase diagram of the URPM 
in the (T, n)-plane. Empty circles: coexistence points ob- 
tained via GCMC simulations (from Ref. [2!|). Full line: fit 
to the coexistence line, assuming a critical exponent /3 = 1. 
Rotated empty square: estimate of the critical point assum- 
ing a critical exponent /3 = 1. Filled circles: spinodal line 
estimated from the condition Snn (& = 0) « 50. Open trian- 
gles: pairing transition points estimated from the condition 
P(m = 1) = 0.3 with r c = 1.4 (as in Ref. [H). Filled trian- 
gles: pairing transition points estimated from the condition 
P(m = 1) « P(m = 2) = 0.5 with r c = 1.0. Stars: CI tran- 
sition points estimated from the vanishing of the electrical 
conductivity (see text for definition). 



are restricted to the symmetric version of the model, 
the URPM. Interactions between polyions are purely 
Coulombic; in particular, contrary to the familiar re- 
stricted primitive model of electrolytes, no hard cores are 
involved. We have reported extensive simulation results 
and theoretical predictions characterizing the pair struc- 
ture, clustering, thermodynamic and dynamical proper- 
ties of the URPM over a wide range of temperatures and 
densities. Concerning the static (structural and ther- 
modynamic) properties, the RPA proves to be a quan- 
titatively reliable approximation at high densities, corre- 
sponding to the fully ionic, conducting regime. At low 
temperature and density, however, the RPA and even the 
HNC theory fail completely, due to the formation of long- 
lived anion/cation pairs, and some larger, mostly neutral 
clusters. 

The low density, low temperature data for pairing, di- 
electric permittivity e and electrical conductivity <r e pro- 
vide strong evidence of a transition between an insulating 
phase, characterized by finite values of e and a vanish- 
ing a e (signaling a vanishing concentration of unpaired 
ions), and a fully ionic, conducting state at higher n and 
T, where a e takes on non-zero values and e is very large. 
Further, indirect evidence for a CI transition is provided 
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by the unusual variation of the polyion mobility with 
temperature at low density. 

At sufficiently low T the CI transition gives way to a 
first order phase separation between a low density insu- 
lating phase and a high density conducting phase. The 
analysis of simulation data for relatively small system 
sizes provides some preliminary evidence for the existence 
of an upper tricritical point terminating the phase coex- 
istence line near the junction with the CI line [29[ . The 
exact nature of the CI transition above the putative tri- 
critical point is not entirely clear. The simulation results 
presented herein and in Ref. [29| point to a continuous 
transition, possibly broadened by finite size effects. The 
tentative phase diagram seem to differ qualitatively from 
that of the RPM, where criticality belongs almost cer- 
tainly to the Ising universality class [H|, Ef| , and there 
is no strong evidence for a CI line above T c . The qual- 
itatively different behaviors of the hard core and pene- 
trable electrolytes are most probably linked to the very 
different nature of the anion/ cation pairs which dominate 
the insulating phase: they are weakly interacting, polar- 
izable entities in the case of the URPM, while they are 
strongly interacting dipolar "molecules" in the RPM (4(| . 
Interestingly, the phase diagram of the URPM is more 
reminiscent of that of the two-dimensional Coulomb gas 
involving logarithmic interactions between ions (oppo- 
sitely charged hard disks). In the low density limit, this 
model is known to undergo an infinite order, Kostcrlitz- 
Thouless (KT) [47| transition between a dielectric phase 
of bound ion pairs and a conducting phase of free ions. 
In the zero density limit the KT transition is character- 
ized by a discontinuous jump of the dielectric permittiv- 
ity from a finite value to infinity as the transition tem- 
perature is approached from below. MC simulations on 
finite systems show that the KT (or CI) transition is 
continuous (rounded), and that the transition tempera- 
ture drops with increasing density; the CI line terminates 
close to the critical point of a first order vapor-liquid co- 
existence curve [ID, |4{|. Within statistical uncertainties, 
the simulations show some evidence of a cusp at the top 
of the vapor- liquid coexistence curve, suggesting a tricrit- 
ical point (rather than the regular critical point observed 
in the three-dimensional equivalent of the Coulomb gas, 
namely the RPM). These findings are confirmed by mean 
field calculations within the chemical representation of a 
mixture of free ions and bound (Bjerrum) pairs (50j . 

In order to confirm a similar scenario in the three- 
dimensional URPM, it will be crucial to investigate fi- 
nite size effects within a full finite size scaling analysis, 
requiring extensive further simulations for several system 
sizes (for a pedagogical review of finite size scaling tech- 
niques, see 5l|). The natural order parameters for such 
an analysis are An, i.e., the difference in density of the 
two coexisting phases, and the fluctuation of the total 
dipole moment per unit volume, which is intimately re- 
lated to the dielectric permittivity (see Eq. (|49l0 : such 
an investigation is under way and will be the main object 
of a forthcoming publication. 



Future developments of our work include an extension 
of the URPM to non-symmetric versions of the UPM, 
which will involve less straightforward aggregation pat- 
terns, and the generalization of the model to oppositely 
charged polyions in the presence of microscopic co- and 
counterions (e.g., from added salt) which will lead to 
screened effective interactions between the polyions as 
opposed to the bare Coulombic interactions considered 
in the present work. The objective of this generaliza- 
tion will be to investigate the influence of screening by 
microions on the CI transition and phase separation. 
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Appendix A: Simulation methods 

Numerical simulations of the URPM have been per- 
formed for systems composed of N = 1000 polyions in 
a cubic cell with periodic boundary conditions using the 
Ewald summation scheme to account for the long range 
of the interactions. In the following, we briefly summa- 
rize the key equations and the parameters employed in 
our simulations. Following a standard procedure [521 ] . we 
introduce screening charge distributions of shape 
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such that the screening charge distribution around a par- 
ticle of species a is —Z a p'(r — r.;) [cf. Eq. fl3J)]. Note 
that, in general, the width a' of the screening distribution 
needs not coincide with the one of actual Gaussian dis- 
tribution. Using the standard Ewald summation scheme, 
the total interaction energy U can then be expressed as 



u = u k + u r - u s 



(A2) 



where Uk and U r are Fourier and real space contributions, 
respectively, and U s is a self-term correcting for the inter- 
action between the screening distribution and its equal 
and opposite compensating distribution. For the URPM, 
the three terms read 
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FIG. 19. (a) Concentration-concentration structure factors 
Scc(k) and (b) number-number structure factors 5*jv]v(fc) at 
n = 0.0035 and different temperatures (as labeled), obtained 
using different values of the Ewald sum precision e: e — 10~ 4 



(full lines), e = 
lines). 



10 3 (dashed lines), and e = 10 (dotted 



where = X)i ^ ex P ' r i) an d we have introduced 
the notation a = \/ a 2 + a' 2 . Note that the familiar ex- 
pression for Uk in the case of purely Coulombic inter- 
action is recovered by substituting a — s- a in Eq. (|A3I) . 
Expressions for the forces between charge distributions 
(needed in MD simulations) are derived in a similar man- 
ner. 

When a' = a the real space contribution is identically 
zero and the calculation of the potential energy and the 
forces can be carried out entirely in Fourier space. Note 
that this is possible only because the effective potential 
Vapir) [see Eq. ©] is bounded. The sum in Eq. (|A3|) is 
carried out for all wave vectors k such that |k| is smaller 
than some cut off value k c . The error in the evaluation 
of Uk can then be roughly estimated as (53j 



e = — 



1 exp 
r 



(A6) 



In our simulations, we used indeed a' = u and adjusted 
k c so as to keep the nominal error e constant at 10 ~ 3 . 
As a consequence, the actual number Nk of wave- vectors 
used for the evaluation of Uk varied as a function of the 
density (as Nk ~ V -1 / 3 ), ranging from 2.5 x 10 2 at n = 
0.35 to 4.2 x 10 4 at n = 0.0035. We found that the 
chosen value of e was sufficient to converge structural and 
dynamic properties below the noise level. To illustrate 
this, in figure [19] we compare results for the structure 
factors SNN(k) and Scc{k) at n = 0.0035 obtained using 
e = 10~ 2 , 10~ 3 , and 10 -4 . We note that the effect of 
reducing e becomes slightly more pronounced at low T. 
We remark that these small discrepancies are, however, 
irrelevant for the purposes of this work. 

At even lower densities (n < 0.001), the number of 



wave-vectors required to keep e constant becomes too 
large for an efficient computation of Eq. (|A3[) . In this 
regime, a simple decimation strategy could be used to 
reduce the number of wave- vectors in each spherical shell 
k±Sk. Alternatively, one could use screening charge dis- 
tributions with width a' ^ a, so as to move part of the 
computational effort into real space. In practice, how- 
ever, we hardly obtained any benefit from using an opti- 
mized a', at least in the range of state parameters inves- 
tigated in the current study. In must be noted, in fact, 
that the range of the pair potential entering Eq. (|A4[) 
is relatively long compared to the complementary error 
function, normally encountered for pure Coulombic inter- 
actions. As a consequence, a large cut off must be em- 
ployed to evaluate the real space contribution U r . There- 
fore, complete splitting of the interaction into real space 
and Fourier space contributions docs not appear partic- 
ularly convenient for the system at hand. 

MD simulations were carried out in the microcanonical 
ensemble (NVE) using the velocity Verlet algorithm with 
a time step 5t = 0.08. Equilibration at the temperature 
of interest was achieved by coupling the system to a mas- 
sively stochastic heat bath [54| . Monte Carlo simulations 
have been performed in the canonical (NVT) and grand- 
canonical (/iVT) ensemble. To accelerate equilibration at 
low density and temperature, we also implemented clus- 
ter moves to displace and insert /delete pairs of ions, as 
described in Ref. In NVT simulations, the maxi- 

mal displacements of standard displacement moves and 
cluster moves were adjusted during equilibration so as to 
keep the acceptance ratio around 30%. In (iVT simu- 
lations, both random and biased insertions/deletions of 
pairs were attempted, together with standard displace- 
ment and cluster moves. Analysis of the acceptance ra- 
tios of Monte Carlo moves as a function of state param- 
eters showed that biased insertions/deletions become fa- 
vorable close to the estimated critical region, as expected. 
Proper equilibration of the system at low density was 
checked by comparing results of different thermal histo- 
ries and the different simulation methods. We remark 
that equilibration may become a serious issue at densi- 
ties and temperatures lower than the ones considered in 
this work. In that case, more efficient simulation meth- 
ods should be used 1371. 



Appendix B: Effective interactions 

In this Appendix a brief derivation is provided of the 
effective interaction between two anion/cation pairs, and 
between a free ion and a pair, in the limit of large separa- 
tions. Consider first the case of two pairs. Let ri and r2 
denote the distances between the two oppositely charged 
polyions in the two pairs, and let R be the vector joining 
the CM's of the pairs. For sufficiently low temperatures 
T, the classical relative vibrations of the two ions in a 
pair are of small amplitude, so that their interaction po- 
tential t>_| (?*) may be replaced by its small r limit (HI 
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i.e., (r) = — /3it 4- 3r 2 with u — Q 2 / '(y/ne 1 'a) and in the Taylor expansion one arrives at 



7 = /3u/6a 2 . The configurational part of the internal 
partition function of an isolated pair is hence 
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The configurational part of the partition function of two 
interacting pairs separated by R is 

q 2 {R) = J dri / dr 2 exp{-^ [v + -(ri) 

+« + -(ra) + «i2(R,r 1 ,r a )]}. (B2) 

If R 3> o one may adopt the point dipole approximation 
for the pair/pair interaction, each carrying a dipole nij = 



Pv 12 (R, ri,r 2 ) 
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Substituting (|B3[) into (|B2|) . and choosing R to be the 
polar axis, q 2 may be cast in the form 



q 2 (R) = 2ttc 
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where (p — <pi — ip 2 . For sufficiently large R the last expo- 
nential in the integrand of (|B4|) may be Taylor expanded. 
Elementary calculations show that all odd power terms 
in the expansion vanish upon carrying out the angular in- 
tegrations. Retaining the zeroth and second order terms 
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The effective pair-pair potential is then given by the 
free energy of the interacting pairs minus the sum of the 
internal free energies of each pair, i.e., 



pv 2 , 2 (R) = Pf 2 {R) 
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valid for R^> a. The "van der Waals" exponent 6 agrees 
with the asymptotic form of the effective potential ex- 
tracted from MC simulations (Fig. [12]). The next contri- 
bution would be 0(1/R 12 ). 

The effective ion/pair potential vi^ 2 (R) may be cal- 
culated along similar lines, starting from the coupling 
between a single ion and a pair of instantaneous dipole 
moment Qr 
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The resulting effective pair potential is found to be 

3?r /cr\ 4 



M r ) = ~tVr 



(B8) 



Note that, contrary to (3v 2 , 2 (R) in Eq. (|B6|) /3vi^{R) de- 
pends on temperature. The same calculation leads to the 
following expression for the reduced polarizability of an 
anion/cation pair 



— ^ = 6\/7re' . 



(B9) 



[1] G. Gouy, J. Phys. 9, 457 (1910). 

[2] D. L. Chapman, Phil. Mag. 25, 475 (1913). 

[3] P. W. Debye and E. Huckel, Phys. Z. 24, 185 (1923). 

[4] N. Bjerrum, Kgl. Danske Vidensk. Selsk. Mat.-Fys. 

Medd. 7, 1 (1926). 
[5] G. Stell, K. C. Wu, and B. Larsen, Phys. Rev. Lett. 37, 

1369 (1976). 

[6] A. Z. Panagiotopoulos, J. Phys.: Condens. Matter 17, 
S3205 (2005). 

[7] J. -P. Hansen and H. Lowen, Ann. Rev. Phys. Chem. 51, 
209 (2000). 

[8] R. van Roij, M. Dijkstra, and J.-P. Hansen, Phys. Rev. 

E 59, 2010 (1999). 
[9] M. E. Leunissen, C. G. Christova, A.-P. Hynninen, C. P. 

Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van 

Roij, and A. van Blaaderen, Nature 437, 235 (2005). 
[10] J. Ryden, M. Ullner, and P. Linse, J. Chem. Phys. 123, 

034909 (2005). 

[11] V. Dahirel and J.-P. Hansen, J. Chem. Phys. 131, 084902 
(2009). 



[12] A.-P. Hynninen, M. E. Leunissen, A. van Blaaderen, and 

M. Dijkstra, Phys. Rev. Lett. 96, 018303 (2006). 
[13] J. B. Caballero, E. G. Noya, and C. Vega, J. Chem. Phys. 

127, 244910 (2007). 
[14] E. Sanz, M. E. Leunissen, A. Fortini, A. van Blaaderen, 

and M. Dijkstra, J. Phys. Chem. B 112, 10861 (2008). 
[15] J.-L. Barrat and J. F. Joanny, Adv. Chem. Phys. 94, 1 

(1995). 

[16] B. Philipp, H. Dautzenberg, K. J. Linov, J. Kotz, and 
W. Davydoff, Prog. Polym. Sci. 14, 91 (1989). 

[17] E. Tsuchida, Pure Appl. Chem. A 31, 1 (1994). 

[18] H. Dautzenberg, Macromolecules 30, 7810 (1997). 

[19] H.-M. Buchhammer, M. Mende, and M. Oelmann, Col- 
loids and Surfaces A: Physicochemical and Engineering 
Aspects 218, 151 (2003). 

[20] S. F. Edwards, Proc. Phys. Soc. London 85, 613 (1965). 

[21] G. H. Fredrickson, The Equilibrium Theory of Inhomo- 
geneous Polymers (Oxford University Press, 2006). 

[22] M. Castelnovo and J. F. Joanny, Eur. Phys. J. E 6, 377 
(2001). 



17 



J. Lee, Y. O. Popov, and G. H. Fredrickson, J. Chem. 
Phys. 128, 224908 (2008). [39 
A. Y. Grosberg, P. G. Khalatur, and A. R. Khoklov, [40 
Mater. Chem. Phys. 3, 709 (1982). 

J. Dautenhahn and C. K. Hall, Macromolecules 27, 5399 [41 
(1994). [42 
P. G. Bolhuis, A. A. Louis, J.-P. Hansen, and E. J. Meijer, [43 
J. Chem. Phys. 114, 4296 (2001). [44 
Y. Levin and M. E. Fisher, Physica A. 225, 164 (1996). 

C. Valeriani, P. J. Camp, J. W. Zwanikken, R. van Roij, [45 
and M. Dijkstra, Soft Matter 6, 2793 (2010). 

D. Coslovich, J.-P. Hansen, and G. Kahl, Soft Matter 7, [46 
1690 (2011). 

M. Baus and J.-P. Hansen, Physics Reports 59, 1 (1980). [47 
F. H. Stillinger and R. Lovett, J. Chem. Phys. 48, 3858 [ 
(1968). 

M. J. Gillan, Mol. Phys. 49, 421 (1983). [49 
D. Ruelle, Statistical Mechanics: Rigorous Results 
(World Scientific Publishing Company, 1999). [50 
J.-P. Hansen and I. R. McDonald, Phys. Rev. A 23, 2041 
(1981). [51 
J.-P. Hansen and I. R. McDonald, Theory of Simple Liq- [52 
uids, Third Edition (Academic Press, 2006). 
J. Kirkwood, J. Chem. Phys. 3, 300 (1935). [53 
C. Valeriani, P. J. Camp, J. W. Zwanikken, R. van Roij, 
and M. Dijkstra, J. Phys.: Condens. Matter 22, 104122 [54 
(2010). 

[38] M. Fisher and S. Fishman, Phys. Rev. Lett. 47, 421 



(1981). 

L. Belloni, J. Chem. Phys. 98, 8080 (1993). 

G. Orkoulas and A. Z. Panagiotopoulos, J. Chem. Phys. 

101, 1452 (1994). 

B. Widom, J. Chem. Phys. 39, 2808 (1963). 

J. Kirkwood, J. Chem. Phys. 7, 911 (1939). 

M. Neumann, Mol. Phys. 50, 841 (1983). 

E. Luijten, M. E. Fisher, and A. Z. Panagiotopoulos, 

Phys. Rev. Lett. 88, 185701 (2002). 

J.-M. Caillol, D. Levesque, and J. -J. Weis, J. Chem. 
Phys. 116, 10794 (2002). 

J. M. Romero-Enrique, L. F. Rull, and A. Z. Pana- 
giotopoulos, Phys. Rev. E 66, 041204 (2002). 
J. Kosterlitz and D. Thouless, J. Phys. C 6, 1181 (1973). 
J.-M. Caillol and D. Levesque, Phys. Rev. B 33, 499 
(1986). 

G. Orkoulas and A. Z. Panagiotopoulos, J. Chem. Phys. 
104, 7206 (1996). 

Y. Levin, X. J. Li, and M. E. Fisher, Phys. Rev. Lett. 
73, 2716 (1994). 

N. B. Wilding, Am. J. Phys. 69, 1147 (2001). 
D. Frenkel and B. Smit, Understanding Molecular Simu- 
lation (Academic Press, 2001). 

W. Smith and T. R. Forester, J. Mol. Graph. 14, 136 
(1996). 

M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Clarendon Press, 1987). 



